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ABSTRACT 

Three-dimensional turbulent jets in crossflow at low to medium jet-to-crossflow 
velocity ratios are computed with a finite- volume numerical procedure which utilizes 
a second-moment closure model to approximate the Reynolds stresses. A multigrid 
method is used to accelerate the convergence rate of the procedure. Comparison of 
the computations to measured data show good qualitative agreement. All trends are 
correctly predicted, though there is some uncertainty on the height of penetration of the 
jet. The evolution of the vorticity field is used to explore the jet-crossflow interaction. 

INTRODUCTION 

Three-dimensional turbulent jets in crossflow have important engineering applications in 
both confined and unconfined environments. Examples of jets issuing into confined crossflow 
include internal cooling of turbine blades, dilution air jets in combustion chambers, jets 
from V/STOL aircraft in transition flight, etc. The examples of turbulent jets issuing into 
unconfined (semi-infinite) crossflow include discharges from cooling towers or chimney stacks 
into the atmosphere, or of sewerage or waste heat into water bodies, film-cooling of turbine 
blades, etc. Several experimental investigations give insight into the various characteristics 
of the jet-crossflow interaction. The earlier studies typified by Kamotani and Greber [1] 
measured mostly global properties such as jet trajectory and spreading rates for a wide 
range of conditions. These results showed that the trajectory based on the temperature 
field often differed from the trajectory based on the velocity field. The same applies to the 
spread rates. There is thus some ambiguity as to what constitutes jet fluid as opposed to 
the jet boundaries and extent of mixing. Flow visualization studies by Foss[2] present an 
insight into the differences in flow regimes between lower and higher jet-to-crossflow velocity 
ratios, R. For the same configurations, measurements of the mean-flow and several higher 
order turbulent statistics are presented by Andreopoulos and Rodi[3] and Andreopoulos[4j. 
These constitute the most detailed set of measured data of jets-in-crossflow existing in the 
open literature. However, in regions with reverse flow, or high turbulence levels the accuracy 
of the data would not be very high. There is therefore still need for “accurate” numerical 
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computations in understanding the flow phenomena. Sykes et al [5] used a numerical method 
to explore the vorticity dynamics of jets in crossflow at high R. Their computations showed 
good qualitative agreement with the measured data [3] at R— 2, except for turbulent kinetic 
energy levels, which were grossly overpredicted in the wake of the jet. 

The present paper examines flow features at moderate R, using a finite- volume based 
multigrid method to compute the mean-flow and a second-moment closure model for the 
Reynolds stresses. 


MATHEMATICAL MODEL 


Mean-Flow 

The time-averaged, three-dimensional, steady state equations governing the turbulent flow 
and heat/concentration transfer may be expressed in cartesian tensor notation as: Continuity 
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with i— 1,2,3 and j— 1,2,3 representing properties in the lateral, vertical and longitudinal 
directions, respectively. y 3 (— y l ,y 2 ,y 3 ) represents the cartesian coordinates; Ui the cartesian 
velocity components; P the pressure and $ the normalized temperature or concentration, p 
is the density, p is the molecular viscosity and A is the thermal or species diffusivity. The 
equations are expanded by using Einstein’s summation rule for repeated indices. -pu-Uj and 
— puppaxe respectively, the Reynolds stresses and heat/concentration fluxes which must be 
determined by a turbulence model before the system of equations can be closed. 

Turbulence Model 

The transport equations for the Reynolds stresses may be expressed as: 

r\ 

Qyj (pUjUiiij ) — p(d{j + Pij + 7r,j — tij) (4) 

where the terms on the right hand side are respectively, diffusive transport, shear production, 
pressure-strain redistribution, and viscous dissipation. Of these, only the shear production 
can be represented exactly. All other terms must be modeled in order to close the system 
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of equations. Dissipation is modeled with the assumption of local isotropy of small scales to 
yield tij = 2/3e<5^, where e is the rate of dissipation of turbulent kinetic energy, k(= \u^u x ). 
Following Launder, Reece and Rodi [ 6 ], hereafter denoted LRR, we utilize a gradient diffusive 
transport model and a quasi-isotropic model for the pressure-strain. The resulting modeled 
equations may be written as: 
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with the empirical coefficients a,/ 3 , 7 ,ci, and c s given by: 

a = 0.7636 — 0.06 /; (3 - 0.1091 + 0.06/; 7 - 0.182; 
C! = 1.5 - 0.50/; c s = 0.22 


Wall-proximity effects are approximated with the function f which takes a value of unity near 
walls decreasing asymptotically to zero in a completely free stream (see Demuren and Rodi 
[7]). The distribution of e is obtained from the solution to the modeled transport equation: 

= + 2 P CtlPii jz ~ Ce2 P% (®) 


Rather than solve the corresponding equations for the turbulent heat/concentration fluxes 
we employ eddy-diffusivity relations as: 
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p t is the eddy viscosity given by : 
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The empirical coefficients in Eq. ( 6 )- ( 8 ) take values: 

c t = 0.15, c M = 0.09, c el = 1.44, c t2 = 1.92, a t = 0.9 
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Boundary Conditions 

Four types of boundary conditions are encountered, namely: inlet, outlet, symmetry and 
walls. Inlet conditions are specified from experimental data. The outlet is an outflow 
boundary requiring no formal specification of conditions. Along symmetry planes the nor- 
mal gradients of all variables are set to zero, and the normal velocity component is also 
zero. The walls are special in that we do not integrate all the way down, rather we use 
the wall- function method [7] to prescribe the values of the dependent variables at the nodes 
immediately adjacent to the walls. The nodes are located such that the flow is in the fully 
turbulent region where the logarithmic law of the wall may be expected to apply. However, 
in the complex flow field under investigation here, this is unlikely to be the case everywhere. 
An improved wall-matching technique, such as that proposed by Degani and Walker[8] would 
probably be more appropriate. The flow which comes out of the jet hole is in general not 
uniform but is distorted to an extent which depends on R. At high R, there is very little 
distortion so uniform exit velocity may not be in much error. However, at low R (e.g. 0.5), 
Andreopoulos and Rodi [3] found that there was hardly any outflow through the front half 
near the stagnation point, and the vertical velocity near the lee side was measured to be 
more than 150 percent of the mean. The best way to treat this case is to solve for the flow 
within the discharge pipe as well, Dernuren and Rodi [9] have found that the imposition of 
a constant total pressure at the discharge plane produced nearly the correct distribution of 
the vertical velocity. This approach is adopted here. 

The boundary conditions for the turbulence quantities are prescribed along near-wall 
nodes based on the assumption of local equilibrium. Hence, we use the equation: 


Pij + 7 tij — 0 


( 9 ) 


which is effectively an algebraic stress model. Further, the assumption of the logarithmic 
law of the wall enables the values of k and e at the near-wall nodes to be related to the 
friction velocity (7, as: 
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where is the von Karman constant with value 0.42, and y w is the normal distance from 
the wall to the respective node. 


SOLUTION PROCEDURE 

Equations (l)-(8) form a closed set which should be solved simultaneously to yield distribu- 
tions of the mean-flow and turbulence quantities. The system of equations are discretized by 
a semi-implicit, finite- volume procedure based on the SIMPLE algorithm of Patankar and 
SpaldingflO]. The original algorithm, as well as proposed improvements [11,12] have poor 
convergence properties, especially for three-dimensional recirculating flows. This deficiency 
is removed in the present work by a multigrid technique which uses SIMPLEC [1 1] merely 
as a relaxation (smoothing) scheme. Shaw et al [1 3j have shown that the SIMPLE algorithm 
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(or extensions thereof) has good smoothing properties. Details of the multigrid procedure 
are provided below. The semi-implicit method involves a linearization which decouples the 
system of equations because the full Jacobian matrices are not evaluated. Although this 
practice produces a scalar system of matrix equations which are computationally cheap to 
solve, the system may become stiff if large source terms have to be treated explicitly. Such 
a situation would arise if the Reynolds stresses in Eq. (2) have large gradients, as one might 
expect in the flows under consideration here, the consequence of that would be massive 
underrelaxation required to ensure stability with the attendant slow convergence rate. The 
stiffness can be reduced considerably by splitting the Reynolds stress into two parts: 
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The first part is treated explicitly. The second part is added to the molecular diffusion 
term and treated implicitly. The modified momentum equation has the form: 
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The resulting system of equations is now very stable and can be solved with standard iteration 
methods. 


Multigrid Procedure 

In the present work the FAS-FMG (full approximation storage-full multigrid) algorithm 
originally developed by Brandt [14] is employed to solve the hydrodynamic equations. The 
present implementation derives from previous works by Demuren [15] and Varika [16], There 
are however significant differences. First, the present method uses a regular grid system with 
no staggering of the velocity nodes relative to the pressure nodes. The expected odd-even 
decoupling problem is overcome by adding a fourth-order artificial dissipation term to the 
pressure gradient. It can be shown that, with a coefficient of unity, this practice is equivalent 
to the so-called “momentum interpolation” method of Peric [17]. However, there is now the 
flexibility to vary the coefficient all the way down to zero, if necessary. The second difference 
is that the system of equations is now solved in a sequential manner as opposed to the 
coupled approach proposed by Vanka. Numerical experiments showed no advantage in using 
the latter in a multigrid procedure, and it can be shown mathematically that it is less stable 
in a single-grid procedure. Further, the decoupled procedure is more easily vectorized. 

The basic steps of the relaxation process are: 

1. Solve the U\ momentum equation using a guessed pressure field. 

2. Then the f/ 2 momentum equation. 
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3. Then the Uz momentum equation, 

4. Compute the velocities on the faces of the control volume, each by linear interpolation 
plus a fourth-order artificial dissipation term, 

5. Then compute the mass source error in each control volume. 

6. Solve a pressure-correction equation to eliminate the mass source errors, and then 
correct the pressures and corresponding velocity components. 

7. If on the finest grid level, solve the equations for k, e, and u^UTj, as the case may be. 

8. Then solve the temperature equation. 

These steps are repeated until convergence on the current finest grid. The multigrid 
process starts with the prolongation of converged solution on a coarse grid to the next finer 
grid. A fixed V-cycle full multigrid algorithm is then utilized with 10 iterations on the 
coarsest grid and 3 on intermediate grids and one iteration on the finest grid. Residuals 
of the momentum equations, as well as the cell-face mass flux equations are restricted onto 
coarser grids for smoothing. Since we employ the full approximation storage variant, which 
is applicable to a non-linear system [14], it is also necessary to restrict the velocities and 
the temperature. Restriction from fine to coarse grid is by 8- point averages. The difference 
between the restricted values and the smoothened solution on the coarse grid is prolongated 
onto the fine grid using trilinear interpolation. An ADI routine is employed to solve the final 
set of algebraic equations for all variables at all grid levels. The underlying algorithm is the 
tri-diagonal matrix algorithm (TDMA) which is known to be recursive, and would thus not 
normally be vectorizable. However, by a change in the data structure we can make all the 
internal loops of the ADI solver vectorizable on the Cray computers. Although we cannot 
remove the recursivity of the algorithm, the change in data structure ensures that all floating 
point operations are in vector form. Typical saving in total CPU time resulting from this 
change alone is of the order of 50% for an average vector length of 35. 

The equations for turbulent quantities and e are solved only on the current finest 
grid. Corresponding operators on coarser grids are calculated using restricted values for 
these quantities. However, the solution process on any fine grid is started with variable 
values prolongated from the converged solution on the immediate coarser grid. 

In order to investigate the characteristics of jets in crossflow at low to medium R, we 
compute two cases with R — 0.5 and 2.0, respectively. These correspond to some of the 
cases for which Andreopoulos and Rodi [3] and Andreopoulos [4] present measured data of 
mean-flow and turbulent quantities. 

RESULTS AND DISCUSSION 

In the computations a non-uniform grid is utilized, with most points concentrated in the 
regions near the jet discharge and immediately downstream. The computational domain ex- 
tends for 6 diameters from the symmetry axis to a free boundary, in the lateral direction, 11 
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diameters from the wall to a free boundary, in the vertical direction, and 22 diameters from 
the crossflow inlet to exit planes, in the streamwise direction. The 3-level multigrid scheme 
has (11,18,23) points on the coarsest level and (38,66,86) points on the finest level, in the 
y l jV 2 ^y 3 directions, respectively. Figure 1 gives a perspective view of the finest grid distribu- 
tion. Typically, 45 fine grid iterations are required for convergence to a normalized residual 
norm of 5xl0“ 4 . With the overheads of the smoothing, restrictions and prolongations, this 
translates to about 120 total work units or 25 minutes of CPU time on the Cray 2 computer. 
It turns out that the time is evenly divided between the mean-flow and the turbulence field 
computations, which suggests that this is a good candidate for a dual-processor machine. 
The inlet plane for the crossflow is located 4D upstream ot the jet hole, and the outlet plane 
18D. The other two boundaries are located far enough away as to have no influence on the 
jet. 

Figure 2 compares the computed streamwise velocity contours, in the symmetry plane, 
for the two cases, and the normalized temperature contours are compared in Fig. 3. Although 
both jets have a reverse flow region on their lee side, their regimes are quite different. The 
jet at low R has a very small core, beyond which it is quickly bent over by the crossflow 
and blends with the wall boundary layer. At higher R the jet has a core which is roughly 
2D long, considerably shorter than the 4-6 D observed in free jets, and is bent over much 
more gradually. An characteristic feature of jet in crossflow is that the temperature and 
velocity contours are noncoincident. The core of the jet, indicated by the temperature field, 
is actually a region with high shear, and the peak streamwise velocity occurs not in the 
jet but in the crossflow, caused by the acceleration of the latter around the jet. Figure 4 
compares the computations of the velocity field at R = 0.5 with measured data [3]. Turbulent 
quantities are compared in Fig. 5. The vertical velocity profiles along the symmetry plane 
show good agreement. It is notable that the flow reversal from upwards to downwards which 
occurs somewhere between ID and 4D is predicted. This is associated with the downwash 
effect of the wake. Streamwise velocity profiles show good agreement between predicted and 
measured values. The lateral velocity in a longitudinal plane passing through the edge of 
the jet is also well predicted. Comparison of turbulent quantities in Fig. 5 indicates that the 
jet may have penetrated deeper into the crossflow than computed, which is surprising if one 
considers the correct prediction of the vertical velocity magnitudes. A plausible explanation 
is that the measurements may represent the averaged effect of large scale coherent structures 
which were observed near the edge of the jet but which the present turbulence model does not 
explicitly account for. In contrast, the comparisons at high R, show in Fig. 6 that we predict 
a higher penetration of the jet. The vertical velocity field is considerably overpredicted. 
Apart from this, the correct trends and magnitudes are predicted. Figure 7 shows similar 
comparisons for the Reynolds stresses. At 4D, the predicted turbulent kinetic energy levels 
are much higher than the measured values. The predictions indicate that there is significant 
dissipation downstream, whereas the experiments show little change in magnitudes so that 
by 6D and 10D predicted and measured levels are similar. That change is attributable 
to convection and diffusion. Sykes et al [5] also overpredicted the turbulence levels. It is 
believed that this is a result of the inadequacy of the dissipation equation in the region with 
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very high shear rates near the jet exit. 

In Figs. 8-10 we examine the evolution of the vorticity field in cross-stream planes, (OD, 
2D and 4D) as the flow progresses downstream. Figure 8 compares, at low and medium R, 
the lateral vorticity components. All along the wall there is a strong boundary layer vorticity 
which has a positive sense. At low R, a horse- shoe vortex is formed around the jet which 
then binds the jet fluid to the boundary layer. The picture is quite different at higher R. 
First, there is the formation of strong negative vorticity in the stagnation region along the 
front of the jet, which is convected above the jet by 2D and is almost completely diffused 
and dissipated by 4D. A horse-shoe vortex forms below it, which is quite strong at 2D, but 
considerably diminishes in strength by 4D, Between the horse- shoe vortex and the boundary- 
layer vortex is a weaker vortex with opposite sense in the wake of the jet. Figure 9 shows 
the vertical vorticity component. In both cases with low- and medium- R we have a similar 
pattern. There is a strong negative vorticity which emanates near the lip of the jet exit hole 
due to the expansion of the jet as it emerges. It is subsequently convected and dispersed 
by the jet-crossflow interaction. This is the dominant feature, though there is a vorticity 
pattern with the opposite sense near the symmetry plane. This pattern is also seen in Fig, 
10 which shows the streamwise vorticity components. It is associated with the downwash 
effect of the wake and leads to a depression of the velocity and temperature contours near 
the symmetry axis. There is a strong negative vorticity produced by the acceleration of the 
crossflow as it bends around the jet. Some of this is then entrained into the jet through 
the wake region. The rest is wrapped around the jet and is dispersed further downstream. 
At low R, there is a much smaller wake region, so this process takes place within a short, 
distance. 

CONCLUDING REMARKS 

The present mathematical procedure predicts correctly the basic features of turbulent jets 
in crossflow at low to medium jet to crossflow velocity ratio. There is some quantitative 
discrepancy between model results and experimental data which is presumed to be partly 
due to model deficiency and inadequacies in the measurement techniques. Further, boundary 
conditions at the jet discharge plane which are difficult to specify because of the crossflow 
interaction have significant effects on the computed flow field. They are rarely ever measured 
in sufficient detail. The differences in the jet- crossflow interaction at low to medium R are 
examined through the evolution of the vorticity field. 
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FIG. 5. Comparison of computations (lines) and experimental 
data (symbols) [3], R— 0.5 turbulence field 
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